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We design fast protocols to separate or recombine two ions in a segmented Paul trap. By inverse 
engineering the time evolution of the trapping potential composed of a harmonic and a quartic term, 
it is possible to perform these processes in a few microseconds without final excitation. These times 
are much shorter than the ones reported so far experimentally. The design is based on dynamical 
invariants and dynamical normal modes. Anharmonicities beyond the harmonic approximation at 
potential minima are taken into account perturbatively. The stability versus an unknown potential 
bias is also studied. 

PACS numbers: 03.67.Lx, 37.10.Ty, 42.50.Dv 


I. INTRODUCTION 

Trapped cold ions provide a leading platform to implement quantum information processing. Separating ion chains 
is in the toolkit of basic operations required. (Merging chains is the corresponding reverse operation so we shall only 
refer to sraaration hereafter.) It has been used to implement two-qubit quantum gates [ij; also to purify entangled 
states [3)0]: teleport material qubits [3|. Moreover, an architecture for processing information scalable to many 

ions could be developed based on shuttling, separating and merging ion crystals in multisegmented traps 0]. 

Ion-chain separation is known to be a difficult operation Q . Experiments have progressed towards lower final exci¬ 
tations and shorter times but much room for improvement still remains Problems identified include anomalous 

heating, so devisiM short-time protocols via shortcuts to adiabaticity (STA) techniques was proposed as a way-out 
worth exploring [l^ . STA methods intend to speed up different adiabatic operations 0[il without inducing final 
excitations. An example of an elementary (fast quasi-adiabatic) STA approach m was already applied for fast chain 
splitting in Q . Here, we design, using a more general and efficient STA approach based on dynamical normal modes 
[l3lIl4|. protocols to effectively separate two equal ions, initially in a common electrostatic linear harmonic trap, into 
a final configuration where each ion is in a different well. The motion is assumed to be effectively one dimensional 
due to tight radial confinement. The external potential for an ion at q is approximated as 

Vext = a{t)q^ + P{t)q^, ( 1 ) 

which is experimentally realizable with state-of-art segmented Paul Traps 0, [l^ . 

Using dynamical normal modes (NM) [l3 . a Hamiltonian will be set which is separable in a harmonic approxi¬ 
mation around potential minima. By means of Lewis-Riesenfeld invariants [l^ we shall design first the approximate 



FIG. 1: (Color online) Scheme of the separation process. At t = 0 (left), both ions are trapped within the same external 
harmonic potential. At final time tf (right), the negative harmonic term, and a quartic term build a double well external 
potential. The aim of the process is to set each of the ions in a different well without any residual excitation. 
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dynamics of an unexcited splitting, taking into account anharmonicities in a perturbative manner, and from that 
inversely find the corresponding protocol, i.e. the a(t) and /3{t) functions. 

The Hamiltonian of the system of two ions of mass m and charge e is, in the laboratory frame. 
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where pi,P 2 are the momentum operators for both ions, qi,q 2 their position operators, and Cc = eo being the 
vacuum permittivity. We use on purpose a c-number notation since we shall also consider classical simulations. The 
context will make clear if c—numbers or g—numbers are required. We suppose that, due to the strong Coulomb 
repulsion, qi > q 2 . By minimizing the potential part of the Hamiltonian V, we find for the equilibrium distance 
between the two ions, d{t), the quintic equation @ 

l3{t)d^{t) + 2a{t)d^{t)-2Cc = 0, (3) 


which will be quite useful for inverse engineering the ion-chain splitting, even without an explicit solution for d(t). 
At t = 0 a single external well is assumed, ,5(0) = 0 and a(0) > 0, whereas in the final double-well configuration 
P{tf) > 0,a(t/) < 0. At some intermediate time ti the potential becomes purely quartic {a{ti) = 0). Our aim is to 
design the functions a{t) and (3{t) so that each of the ions ends up in a different external well as shown in Fig. [U in 
times as short as possible, and without any final excitation. 


II. DYNAMICAL NORMAL MODES 


To define dynamical NM coordinates, we calculate first at equilibrium (the point in configuration space 

where the potential is a minimum, dV/dqi = dVfdq 2 = 0) the matrix Vij = — Y I . The equilibrium positions 
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The eigenvalues are 


A_ = —(2a-k3^dO, 
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which define the NM frequencies as = y^Ai corresponding to center-of-mass (—) and relative (stretch) motions 
(-k). These relations, with Eq. ([3|) written as 


Pit) = 
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allow us to write a{t) and d{t) as functions of the NM frequencies, 
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a{t) = —m (3H^ — 5H?_) 
d{t) = 


4Cc 


Substituting these expressions into Eq. 
The normalized eigenvectors are 


m (n^ — fl?.) 

/3(t) may also be written in terms of NM frequencies. 
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which we denote as v± = . The (mass weighted) dynamical NM coordinates are defined in terms of the laboratory 

coordinates as 


q± = a±^/m{ql - + b±y/m{q2 - qi°^). 

The unitary transformation of coordinates is 


U = 


j (iq+dq_dgidg 2 |q+, q 
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where (q+, q_|gi, 92 ) = — gi(q+, q-)]i 5 [g 2 — ( 72 (q+,q-)]- The Hamiltonian in the dynamical equation for |'0') = 

where |'0) is the lab-frame time-dependent wave function evolving with H, is given by 
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plus qubic and higher order terms in the potential that we neglect by now (they will be considered in Sec. IV 
below). Similarly to [ij, Hli apply a further unitary transformation U = e“®v^<^q+/('/ 2 fi)^ write down an 
effective Hamiltonian for |'0") = with the form of two independent harmonic oscillators in NM space, H" = 

UH'U"^ - iWKptW), 


P- ^ ^o2„2 
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These oscillators have dynamical invariants of the form 

1 , 


a = r[/)±(p± - i±) - p±(q± - x±)p 

+ in? 


‘0± 


V p± 


(14) 


where the auxiliary functions p± and x+ satisfy 
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with Hgi = Hj-(O), and, due to symmetry, X- = 0. 

The physical meaning of the auxiliary functions may be grasped from the solutions of the time-dependent 
Schrodinger equations (for each NM Hamiltonian H'^ in Eq. (1131) 1 proportional to the invariant eigenvectors [2l|. 
They form a complete basis for the space spanned by each Hamiltonian H'^ and take the form 


(q±|V^„±(t)) = e 


_ i ^nicr±) 
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where a± = a-nd $„(cr±) are the eigenfunctions of the static harmonic oscillator at time t = 0. Thus p± are 

scaling factors proportional to the state “width” in NM coordinates, whereas the x± are the dynamical-mode centers 
in the space of NM coordinates. Within the harmonic approximation there are dynamical states of the factorized 
form |'0"(t)) = |V’4.(t))|'0” (t)) for the ion chain dynamics, where the NM wave functions \ilj±{t)) evolve independently 
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FIG. 2: (Color online) Evolution of a) b) /3(t); and c) d(t). In d) the NM frequencies, solid line for the and dashed line 
for the ‘+’ are depicted. Two ^Be"*" ions were separated in the simulation, with ujo/{2tv) — 2 MHz, tf = 5.2 /rs, ao = ma;o/2, 
and d(0) = 5.80 fj,m. 


with H'^. They may be written as combinations of the form |'0(^(t)) = with constant amplitudes 

Cn±- The average energies of the n-th basis states for the two NM are , 
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III. DESIGN OF THE CONTROL PARAMETERS 

Once the Hamiltonian and Lewis-Riesenfeld invariants are defined, we proceed to apply the invariant based inverse 
engineering technique and design shortcuts to adiabaticity. The results for the simple harmonic oscillator in serve 
as a reference but have to be extended here since the two modes are not really independent from the perspective of 
the inverse problem. This is because a unique protocol, i.e., a single set of a{t) and /3(t) functions has to be designed. 

We first set the initial and target values for the control parameters a{t) and At time t = 0, the external 

trap -for a single ion- is purely harmonic, with (angular) frequency ujq- From Eq. ([5|), we find that II_(0) = luq and 

f2+(0) = a/Swq. The equilibrium distance is d{0) = ^ For the final time, we set a tenfold expansion of the 

equilibrium distance, d{tf) = 10(i(0), and = Wq. This also implies H+(ty) = \/1.002 wq Ri H_(t/), i.e., the hnal 

frequencies of both NM are essentially equal, the Coulomb interaction is negligible, and the ions can be considered to 
oscillate in independent traps. 

The inverse problem is somewhat similar to the expansion of a trap with two equal ions in Ref. 0 , but complicated 
by the richer structure of the external potential. The Hamiltonian (I13II and the invariant (1141) must commute at both 
boundary times [E[{tb), I(tb)] = 0, 


tb = 0,tf, (19) 

to drive initial levels into final levels via a one-to-one mapping. This is achieved by applying appropriate boundary 
conditions (BC) to the auxiliary functions p±, x± and their derivatives, 


p±(0) 

= 1, P±(i/)=7±, 

(20) 

p±{tb) 

= P±{tb) = 0, 

(21) 

X+{tb) 

= x+{tb) = x+{tb) = 0, 

(22) 


where 7 ± = Let us recall that a;_ = 0 for all times so this parameter does not have to be considered further. 

Inserting the BC for x+{tb) and x+{tb) in Eq. (ITbll we find that d{tb) = 0. Additionally, d{tb) = 0 is to be imposed 
so that U{tb) = 1. According to Eq. ([5|), d(4) = d{tb) = 0 by imposing ^±{tb) = ii±(tb) = 0. With d{tb) = d{tb) = 0 
the Hamiltonians and wave functions coincide at the boundaries, H'{tb) = H"{tb), |V'^(4)) = 14’”(db)), which simplifies 
the calculation of the excitation energy. 
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From Eq. m, the NM frequencies may be written as 


fl± 
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(23) 


Thus the BC Ct±{ti,) = = 0 are satisfied by imposing on the auxiliary functions the additional BC 

'P±{tb) = "P±{tb) = 0. (24) 

We may now design ansatzes for the auxiliary functions p± that satisfy the ten BC in Eqs. (I20I21I24L plus the BC 
for x+{tb) and x+{tb) in Eq. (|22ll (since d{tb) = 0, x+{tb) = 0 is then automatically satisfied, see Eq. (|T6)) '). Finally, 
from the NM frequencies given by Eq. (I23|) we can inverse engineer the control parameters a(t) and f3(t) from Eqs. 
(I7I8I) and (0. 

A simple choice for p-{t) is a polynomial ansatz of 9th order p_ = where s = t/tf. Substituting this 

form in the ten BC in Eqs. (I20I21I24L we finally get 

p- = 126(7_ - l)s® - 420(7_ - 1)5^ + 540(7_ - l)s^ 

- 315(7_ - l)s® + 70(7_ - 1)5^ + 1. (25) 

For p+ we will use an 11-th order polynomial to satisfy as well x+{tb) = x+{tb) = 0. The parameters 

ao -9 are fixed so that the 10 BC for p^ are fulfilled (see the Appendix), whereas aio, an are left free, and will be 
numerically determined by a shooting program (‘fminsearch’ in MATLAB which uses the Nelder-Mead simplex 
method for optimization), so that the remaining BC for x+{tb) and x+(tb) are also satisfied. Specifically, for each pair 
{aio,aii}, and d{t) are determined from Eqs. (fTKl) and ([8]), to solve Eq. (fTHjl for x+{t) with initial conditions 

x+(0) = i+(0) = 0. The free constants are changed until x+{tf) = 0 and x+{tf) = 0 are satisfied. Numerically a 
convenient way to find the solution is to minimize the energy in Eq. (11811 . 

Eig. (a) and (b) depict the control parameters a{t) and /3(t) found with this method, using Eqs. ([7]) and ([6|), 
for some value oft/ and ojq, see the caption, while Eig. [2] (c) represents the equilibrium distance between ions as a 
function of time ([5]), and Fig. [l](d) the NM frequencies. In Fig. [3](a) the excitation energy is shown, versus final time, 
for optimized parameters given in Fig. (b). The initial state is the ground state of the two ions. It is calculated 
with an “imaginary time evolution”. This is a variation of any numerical time evolution method (here we used a 
variation of the split-operator method), where instead of the real time one defines imaginary time in the evolution 
operator. By letting an ansatz wave function evolve for the static initial potential, it will eventually converge to the 
ground state of the system. The excitation energy is = E{tf) — Eo{tf), where E{tf) is the final energy, calculated 
in the lab frame, and E^itf) is the final ground-state energy. The wave function evolution is calculated using the 
“Split-Operator Method” with the full Hamiltonian ([2|). If the harmonic approximation were exact, there would not 
be any excitation with this STA method, E{tf) = E'^j^itf) + EQ_{tf) = Eo(t/), see Eq. (fT^ . The actual result is 
perturbed by the anharmonicities and NM couplings. The final ground state is also calculated with an “imaginary 
time evolution”. The corresponding final ground state energy is essentially two times a harmonic oscillator ground 
state energy plus the (negligible) Coulomb repulsion at distance d{tf). Eor the final times of all the examples, as it 
was noted in previous works classical simulations (solving Hamilton’s equations from the equilibrium 

configuration instead of Schrodinger’s equation) give indistinguishable results in the scale of Eigure[3] (a). 

The excitation energy in Figure [31(a) (solid line) increases at short times since the harmonic NM approximation fails 
[r3.[l^. However, it goes down rapidly below one excitation quantum at times which are still rather small compared 
to experimental values used so far [3, S] ■ In the following section we shall apply a perturbative technique to minimize 
the excitation further. 


IV. BEYOND THE HARMONIC APPROXIMATION 

An improvement of the protocol is to consider the perturbation of the higher order terms neglected in the Hamil¬ 
tonian 031). These “anharmonicities” 0 are cubic and higher order terms in the Taylor expansion of the Coulomb 
term Cc/{qi - 92 ), 
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FIG. 3: (Color online) a) Final excitation energy vs. final time using the inverse engineering design of Sec. Ill (solid blue), 
and the design that takes into account anharmonicities in Sec. IV (dashed red), b) Values of the free parameters aio (solid 
blue) and an (dashed red) that minimize the excitation energy for the 11th order polynomial (EH). c) Parameters cio (solid 
blue), cii (dashed red) and ci 2 (dash-dotted green) that minimize the excitation energy for the 12th order polynomial (IA2II . 
Two ^Be"*" ions where split, with ujq/{2'k) = 2 MHz. 


Wo (MHz) 

Pmax (10“®N/m®) 

tcrit (M®) 

3 

44.2 

2.9 

2 

11.4 

4.4 

1.2 

2.082 

7.4 

0.8 

0.539 

11.2 


TABLE I: Maximum values of /3, and critical times (final times at which excitations below 0.1 quanta are reached) for different 
values of ujq. The calculations were performed with the 11-th order polynomial for p+. 


In NM coordinates the terms take the simple form 


( 27 ) 

which may be regarded as a perturbation to be added to iL" in Eq. ([13]). (The perturbation does not couple the 
center-of-mass and relative subspaces.) To first order, the excess energy due to these perturbative terms at final time 
is given by 

= (28) 

where the \'4>n+) the unperturbed states in Eq. (HID. Inverse engineering the splitting process may now be carried 
out by considering a 12th order polynomial for (see (|A2p ). with three free parameters so as to fix the BC for 
a;+ and also minimize the excitation energy. In practice we use MATLAB’s ‘fminsearch’ function for the shooting to 

minimize Eo+{tf) +SEq^ as no significant improvement occurs by including higher order terms. Fig. [3] (a) compares 
the performance of such a protocol with the simpler one with the llth-order polynomial (lAll) . Fig. |3| (c) gives the 
values of optimized parameters at different final times. 


V. DISCUSSION 

A large quartic potential is desirable to control the excitations produced at the point where the harmonic term 
changes its sign |10l |. At this point, the harmonic potential switches from confining to repulsive, which reduces the 
control of the system and potentially increases diabaticities and heating. In the inverse approach proposed here there 
is no special design of the protocol at this point, but the method naturally seeks high quartic confinements there. 
In Fig. [2] (b) /3 reaches its maximum value right at the time where a changes sign (see Fig. [2] (a)). However, the 
maximum value that /3 can reach will typically be limited in a Paul trap i. In Table m we summarize the different 
maximal values of /3 and critical times (final times at which excitations below O.I quanta are reached) for different 
values of loq using the 11-th order polynomial (lAll) for p_)_. The maximum P decreases with tf, such that the shortest 
possible t/ at a given maximum tolerable excitation energy is limited by the achievable /3. The trap used in Ref. Q 
yields a maximum P of about 10“"^ N/m^, at ±10 V steering range. In a recent experiment reported in where 
although the purpose was not ion separation a double well potential was produced, the value used was /3 ~ 5 x 10“^ 
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FIG. 4: (Color online) Excitation energy vs. different tilt values of the external potential in terms of the energy difference 
between both wells (upper axis) and values of the A parameter (lower axis), when using the 11-th order polynomial in the 
evolution. Same parameters as in Fig. [2] 



FIG. 5: (Color online) Excitation energy vs. final time comparing the 11-th order polynomial (solid blue) and a non optimized 
trajectory experimentally used in (dashed red) in the evolution. Same parameters as in Fig. [21 


N/m^. The numbers reported in the Table are thus within reach, as the /3 coefficients scale with the inverse 4th 
power of the overall trap dimension, and technological improvements on arbitrary waveform generators may allow for 
operation at an increased voltage range. 

Another potential limitation the method could encounter in the laboratory is due to biases (a linear slope) in the 
trapping potential, Vext = ctq^ + j3q^ + Xq, with A constant and unknown [l3 |. Fig. |4] represents the excitation energy 
versus the energy difference between the two final minima of the external potential, AE (also vs. A). To calculate the 
results, a(t) and /3(t) are designed as if A = 0, but the dynamics is carried out with a non-zero A, in particular the 
initial state is the actual ground state, including the perturbation. Note that AE should be more than a thousand 
vibrational quanta to excite the final energy by one quantum. In Ref. Q an energy increase of ten phonons at about 
150 zN and 80 ^s separation time was reported, so the STA ramps definitely improve the robustness against bias. 

Further experimental limitations may be due to random fluctuations in the potential parameters, or higher order 
terms in the external potential. We leave these important issues for a separate study but note that the structure of 
the STA techniques used here is well adapted to deal with noise or perturbations [231 - 1^ . 

Finally, we compare in Fig. [5] the performance of the protocols based on the polynomials (1251) and (lAll) with a 
simple non-optimized protocol based on those experimentally used in [^. There, the equilibrium distance d is first 
designed as d{t) = d{0) -|- [d{tf) — c?(0)]s^ sin(s7r/2), where s = t/tf. From the family of possible potential ramps 
consistent with this function, we chose a polynomial that drives a from a(0) = oq to a(t/) = —ao/2 (as in Fig. [^l 
and whose hrst derivatives are 0 at both boundary times. /3 is given by Eq. For the times analysed in Fig. 
El the method based on Eqs. (EH) and dHI) clearly outperforms the non-optimized ramp. To get excitations below 




the single motional excitation quantum with the non-optimized protocol, final times as long as t/ ~ 80/xs would be 
needed, which is in line with current experiments. 

We conclude that the method presented here, could bring a clear improvement with respect to the best results 
experimentally reported so far [3, S]- The parameters required are realistic in current trapped ions laboratories. The 
simulations show that, under ideal conditions, the separation of two ions could be performed in a few oscillation 
periods, at times similar to those required for other operations as transport 0 or expansions 0 , also studied with 
STA. 
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Appendix A: Ansatz for p+ 


The ansatz for that satisfies the BC / 9 +( 0 ) = 1, p+{tf) = 7 +, p+{tb) = p+{tb) = P + {tb) 
free parameters takes the form 


= 1 — (126 — 1267 + -I- oio + 5aii)s® 
-k (420 - 42 O 7 + -f 5aio -f 24aii)s® 

— (540 — 54 O 7 + -I- lOoio + 45aii)s^ 
-f (315 - 3157 +-flOaio-f40aii)s® 

- (70- 707 +-f 5aio + 15aii)s® 

-I- aios^° + 


To minimize the perturbation energy in Eq. (I28p . three free parameters are introduced, 

p+ = 1 — (126 — 1267+-I-cio + 5 cii-I- 15 ci 2 )s^ 

-f (420 — 42 O 7 - 1 - + 5cio + 24cii -I- 70ci2)s^ 

- (540 - 54 O 7 + -f lOcio + 45cii -k 126 ci 2 )s^ 

-k (315 - 3157 + -k lOcio + 40cii -k 105 ci2)s® 

— (70 — 7 O 7 + -k 5cio + 15cii -k 35 ci2)s® 

-k CioS^° -k -k Ci2S^^. 


P +{tb) =0 with two 
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